function [AUC,minMC,VAL]=roc(obs,pred,plt,step)

% ROC Response Operating Characteristic curves
% Syntax: [AUC,minMC,VAL]=roc(obs,pred,{plt},{step})
% Description:
%    Creates a ROC curve plot of sensitivity vs 1-specificity
%    which is equivalent of p(true positives) vs p(false negatives)
%
%    Input Arguments:
%
%    obs  -  observed outcome
%    pred -  predicted outcome
%    plt:
%        1 - regular ROC curve plot
%        2 - ROC curve with threshold values
%        3 - Show also false discovery rate
%            any other option will cause the graphic
%            display not to be produced. Default value is 0.
%    step  - this argument controls the stepwiese changes in threshold
%            value. The default value is to use every single outcome value
%            predicted so the only reason to change this argument is if
%            the number of such values is so large that it is affecting the
%            performance of your application. Note this argument can by a
%            single number or a vector, depending if fixed steps are
%            desired or a more flexible specification fo all steps is
%            needed.
%
%
%    Output arguments
%
%    AUC - area under the curve
%    VAL - values defining the ROC curve
%
% Jonas Almeida, May 2006

if nargin<3;plt=0;end % By default do not plot
if(nargin<4 )
    step=unique(pred)';
elseif length(step)==1 % the a fixed step size is being requested
    step=[min(pred):step:max(pred)];
end

obs_pos=sum(obs==1); % Observed positives = 100% true positives
positivos=pred(find(obs==1));npos=length(positivos);
negativos=pred(find(obs==0));nneg=length(negativos);


%Graphic display of ROC curve
switch plt
    case 0
        VAL=[];
        for tr=step    % tr= threshold
            pred_pos = (pred>=tr);pred_neg = (pred<=tr);
            mc=(pred_pos~=obs); % miss-classification: false positives and false negatives
            VAL=[VAL;[tr,sum(negativos>=tr),sum(positivos>=tr),sum(mc),sum(pred_pos.*(1-obs))/sum(pred_pos),sum(pred_neg.*(obs))/sum(pred_neg)]];%sum(mc.*obs)/sum(pred>=tr),1]];%sum(mc.*(1-obs))/sum(pred<=tr)]];
        end
        VAL(:,2)=VAL(:,2)./nneg; %true negatives
        VAL(:,3)=VAL(:,3)./npos; %true positives
        VAL(:,4)=VAL(:,4)./(nneg+npos); % miss-classification: false positive and negative discovery rate
        % inverting order
        VAL=VAL(end:-1:1,:);

        % Area under the curve:
        arestas=[VAL(1:end-1,2),VAL(2:end,2),VAL(1:end-1,3),VAL(2:end,3)]; % edges
        AUC=sum((arestas(:,2)-arestas(:,1)).*(arestas(:,3)+(arestas(:,4)-arestas(:,3))./2)); %area under the ROC curve
        %FDR
        [MCmin,Ind_MCmin]=min(VAL(:,4));
        minMC=VAL(Ind_MCmin,:);
    case 1
        [AUC,minMC,VAL]=roc(obs,pred,0);
        [MCmin,Ind_MCmin]=min(VAL(:,4));
        h=area(VAL(:,2),VAL(:,3));
        set(h,'FaceColor',[0.75 0.75 0.75])
        hold on
        plot([0;1],[0;1],':k');
        ylabel('sensitivity (true positives)')
        xlabel('1-specificity (true negatives)')
        text(0.6,0.4,[' Area = ',num2str(AUC)],'FontSize',14)
        gridover
        hold off

    case 2 % plot threshold value too        
        [AUC,minMC,VAL]=roc(obs,pred,0);
        [MCmin,Ind_MCmin]=min(VAL(:,4));
        h=area(VAL(:,2),VAL(:,3));hold on
        set(h,'FaceColor',[0.75 0.75 0.75])
        set(h,'LineWidth',2)
        set(gcf,'Color','w') %white background
        %plot(VAL(:,2),VAL(:,3),'-','LineWidth',2)
        hy=ylabel('   Sensitivity (true positives) [area],\newline miss-classification [blue], FDR [red]');
        set(hy,'FontName','Verdana','HorizontalAlignment','center')
        hx=xlabel('1-specificity (true negatives)');
        set(hx,'FontName','Verdana')
        text(1,1,[' AUC = ',num2str(AUC),'  '],'FontSize',10,'Color','k','HorizontalAlignment','right','VerticalAlignment','bottom','FontName','Verdana','FontAngle','italic','FontWeight','bold')
        set(get(h,'Parent'),'XMinorGrid','on','YMinorGrid','on')
        hold on
        gridover
        %plot([0;1],[0;1],'--k','LineWidth',2,'Color',[0.2 0.2 0.2]);
        hold off
        ax=[0 1 0 1];
        % Plot FDR
        hold on
        plot(VAL(:,2),VAL(:,5),'LineWidth',2,'Color','r');
        % and the miss-classification
        hfdr=plot(VAL(:,2),VAL(:,4),'LineWidth',2);
        %plot(VAL(Ind_MCmin,2),MCmin,'o','MarkerFaceColor','b'); % put circle on minimum FDR
        %angulo=0;%20;
        %htext=text(VAL(Ind_MCmin,2),MCmin,sprintf('%.2f   ',MCmin),'HorizontalAlignment','right','VerticalAlignment','top','Color','b','FontName','Arial','FontWeight','bold','FontSize',12,'FontAngle','italic','Rotation',angulo);
        %if VAL(Ind_MCmin,2)<0.2
        %    set(htext,'HorizontalAlignment','left','String',['  MCmin=',sprintf('%.2f',MCmin)],'Rotation',-angulo)
        %end
        %htext=text(VAL(Ind_MCmin,2),VAL(Ind_MCmin,3),[sprintf('(%.2f,%.2f) ',VAL(Ind_MCmin,2:3)),' '],'HorizontalAlignment','right','VerticalAlignment','bottom','Color','b','FontName','Arial','FontWeight','bold','FontSize',12,'FontAngle','italic','Rotation',-angulo,'Color','k');
        %if VAL(Ind_MCmin,2)<0.2
        %    set(htext,'HorizontalAlignment','left','String',[sprintf('(%.2f,%.2f)',VAL(Ind_MCmin,2:3))],'Rotation',-angulo)
        %end
        %plot(VAL(:,2),VAL(:,5),'LineWidth',1,'Color','r');
        %plot(VAL(:,2),VAL(:,6),'LineWidth',1);
        hold off
        % Plot Threshold
        hthr=axes;
        axes(hthr);
        plot(VAL(:,2),VAL(:,1),'LineWidth',2,'Color',[0 0.5 0]);
        set(hthr,'YAxisLocation','right','Color','none')
        ylabel('Segmentation threshold','FontName','Verdana','Rotation',-90,'Color',[0 0.5 0],'VerticalAlignment','bottom');
        hold on
        plot(VAL(Ind_MCmin,2),VAL(Ind_MCmin,1),'o','MarkerFaceColor',[0 .5 0],'MarkerEdgeColor',[0 .5 0],'MarkerSize',5); % put circle on minimum FDR
        text(minMC(2),minMC(1),' \rightarrow','HorizontalAlignment','left','Color',[0 0.5 0],'FontWeight','bold');
        text(minMC(2),minMC(1),[num2str(minMC(1)),'   '],'HorizontalAlignment','right','Color',[0 0.5 0],'FontName','Verdana','FontSize',9,'FontAngle','italic','FontWeight','bold');
        %htext=text(VAL(Ind_MCmin,2),VAL(Ind_MCmin,1),[sprintf('   %.3f',VAL(Ind_MCmin,1))],'HorizontalAlignment','left','VerticalAlignment','bottom','Color',[0 .5 0],'FontName','Arial','FontWeight','bold','FontSize',12,'FontAngle','italic','Rotation',angulo);
        hold off
        % Line for minimum of miss-classification
        %axes (h);
        H=get(h,'Parent');
        axes(H);hold on;
        plot([minMC(2) minMC(2)],[0,1],':','LineWidth',2);
        %Plot values on that line
        plot(minMC(2),minMC(3),'ko','MarkerFaceColor','k','MarkerSize',5); % put circle on minimum miss-classification
        plot(minMC(2),minMC(5),'ro','MarkerFaceColor','r','MarkerSize',5); % put circle on corresponding FDR
        plot(minMC(2),minMC(4),'bo','MarkerFaceColor','b','MarkerSize',5); % put circle on corresponding FDR
        %Place arrows indicating axis
        text(minMC(2),minMC(3),'\leftarrow ','HorizontalAlignment','right','Color','k','FontWeight','bold','FontSize',14);
        text(minMC(2),minMC(3),['   ',num2str(minMC(3))],'HorizontalAlignment','left','Color','k','FontName','Verdana','FontSize',9,'FontAngle','italic','FontWeight','bold');
        text(minMC(2),minMC(5),'\leftarrow ','HorizontalAlignment','right','Color','r','FontWeight','bold','FontSize',14);
        text(minMC(2),minMC(5),['   ',num2str(minMC(5))],'HorizontalAlignment','left','Color','r','FontName','Verdana','FontSize',9,'FontAngle','italic','FontWeight','bold');
        text(minMC(2),minMC(4),'\leftarrow ','HorizontalAlignment','right','Color','b','FontWeight','bold','FontSize',14);        
        text(minMC(2),minMC(4),['   ',num2str(minMC(4))],'HorizontalAlignment','left','Color','b','FontName','Verdana','FontSize',9,'FontAngle','italic','FontWeight','bold');
        %Write values of intercept
        text(minMC(2),1,num2str(minMC(2)),'VerticalAlignment','bottom','HorizontalAlignment','center','Color','b','FontName','Verdana','FontSize',9,'FontAngle','italic','FontWeight','bold');
        
        %return second axes to the top:
        axes(hthr);
        %Legend
        % FDR plot
        f=figure;
        set(2,'Color','w');
        plot(VAL(:,1),VAL(:,5),'r','LineWidth',2);hold on
        plot(VAL(:,1),VAL(:,6),'LineWidth',2,'Color',[0.5 0.5 0.5]);
        ax=axis;
        plot([ax(1),ax(2)],[minMC(6),minMC(6)],':','LineWidth',2,'Color',[0.5 0.5 0.5]);
        plot([minMC(1),minMC(1)],[ax(3),ax(4)],':','LineWidth',2,'Color',[0 0.5 0]);
        plot([ax(1),ax(2)],[minMC(5),minMC(5)],':r','LineWidth',2);
        plot(minMC(1),minMC(5),'ko','MarkerFaceColor','r','MarkerSize',5);
        plot(minMC(1),minMC(6),'ko','MarkerFaceColor',[0.5 0.5 0.5],'MarkerSize',5);
        text(ax(1),minMC(5),[num2str(minMC(5)),'  '],'HorizontalAlignment','right','Color','r','FontName','Verdana','FontSize',9,'FontWeight','bold');
        text(ax(2),minMC(6),['  ',num2str(minMC(6))],'HorizontalAlignment','left','Color',[0.5 0.5 0.5],'FontName','Verdana','FontSize',9,'FontWeight','bold');
        text(minMC(1),ax(4),[' ',num2str(minMC(1))],'VerticalAlignment','top','HorizontalAlignment','left','Color',[0 0.5 0],'FontName','Verdana','FontSize',9,'FontWeight','bold');
        axis(ax);
        legend('FDR','FNR',-1)
        grid on
        xlabel('Threshold value','Color',[0 0.5 0])
        title('Positive (FDR) and negative (FNR) discovery rates')
    case 3
        [AUC,minMC,VAL]=roc(obs,pred,0);
        [MCmin,Ind_MCmin]=min(VAL(:,4));
        h=area(VAL(:,2),VAL(:,3));hold on
        set(h,'FaceColor',[0.75 0.75 0.75])
        set(h,'LineWidth',2)
        set(gcf,'Color','w') %white background
        %plot(VAL(:,2),VAL(:,3),'-','LineWidth',2)
        hy=ylabel('   Sensitivity (true positives) [area],\newline miss-classification [blue], FDR [red]');
        set(hy,'FontName','Verdana','HorizontalAlignment','center')
        hx=xlabel('1-specificity (true negatives)');
        set(hx,'FontName','Verdana')
        text(1,1,[' AUC = ',num2str(AUC),'  '],'FontSize',10,'Color','k','HorizontalAlignment','right','VerticalAlignment','bottom','FontName','Verdana','FontAngle','italic','FontWeight','bold')
        set(get(h,'Parent'),'XMinorGrid','on','YMinorGrid','on')
        hold on
        gridover
        %plot([0;1],[0;1],'--k','LineWidth',2,'Color',[0.2 0.2 0.2]);
        hold off
        ax=[0 1 0 1];
        % Plot FDR
        hold on
        plot(VAL(:,2),VAL(:,5),'LineWidth',2,'Color','r');
        % and the miss-classification
        hfdr=plot(VAL(:,2),VAL(:,4),'LineWidth',2);
        %plot(VAL(Ind_MCmin,2),MCmin,'o','MarkerFaceColor','b'); % put circle on minimum FDR
        %angulo=0;%20;
        %htext=text(VAL(Ind_MCmin,2),MCmin,sprintf('%.2f   ',MCmin),'HorizontalAlignment','right','VerticalAlignment','top','Color','b','FontName','Arial','FontWeight','bold','FontSize',12,'FontAngle','italic','Rotation',angulo);
        %if VAL(Ind_MCmin,2)<0.2
        %    set(htext,'HorizontalAlignment','left','String',['  MCmin=',sprintf('%.2f',MCmin)],'Rotation',-angulo)
        %end
        %htext=text(VAL(Ind_MCmin,2),VAL(Ind_MCmin,3),[sprintf('(%.2f,%.2f) ',VAL(Ind_MCmin,2:3)),' '],'HorizontalAlignment','right','VerticalAlignment','bottom','Color','b','FontName','Arial','FontWeight','bold','FontSize',12,'FontAngle','italic','Rotation',-angulo,'Color','k');
        %if VAL(Ind_MCmin,2)<0.2
        %    set(htext,'HorizontalAlignment','left','String',[sprintf('(%.2f,%.2f)',VAL(Ind_MCmin,2:3))],'Rotation',-angulo)
        %end
        %plot(VAL(:,2),VAL(:,5),'LineWidth',1,'Color','r');
        %plot(VAL(:,2),VAL(:,6),'LineWidth',1);
        hold off
        % Plot Threshold
        hthr=axes;
        axes(hthr);
        plot(VAL(:,2),VAL(:,1),'LineWidth',2,'Color',[0 0.5 0]);
        set(hthr,'YAxisLocation','right','Color','none')
        ylabel('Segmentation threshold','FontName','Verdana','Rotation',-90,'Color',[0 0.5 0],'VerticalAlignment','bottom');
        hold on
        plot(VAL(Ind_MCmin,2),VAL(Ind_MCmin,1),'o','MarkerFaceColor',[0 .5 0],'MarkerEdgeColor',[0 .5 0],'MarkerSize',5); % put circle on minimum FDR
        text(minMC(2),minMC(1),' \rightarrow','HorizontalAlignment','left','Color',[0 0.5 0],'FontWeight','bold');
        text(minMC(2),minMC(1),[num2str(minMC(1)),'   '],'HorizontalAlignment','right','Color',[0 0.5 0],'FontName','Verdana','FontSize',9,'FontAngle','italic','FontWeight','bold');
        %htext=text(VAL(Ind_MCmin,2),VAL(Ind_MCmin,1),[sprintf('   %.3f',VAL(Ind_MCmin,1))],'HorizontalAlignment','left','VerticalAlignment','bottom','Color',[0 .5 0],'FontName','Arial','FontWeight','bold','FontSize',12,'FontAngle','italic','Rotation',angulo);
        hold off
        % Line for minimum of miss-classification
        %axes (h);
        H=get(h,'Parent');
        axes(H);hold on;
        plot([minMC(2) minMC(2)],[0,1],':','LineWidth',2);
        %Plot values on that line
        plot(minMC(2),minMC(3),'ko','MarkerFaceColor','k','MarkerSize',5); % put circle on minimum miss-classification
        plot(minMC(2),minMC(5),'ro','MarkerFaceColor','r','MarkerSize',5); % put circle on corresponding FDR
        plot(minMC(2),minMC(4),'bo','MarkerFaceColor','b','MarkerSize',5); % put circle on corresponding FDR
        %Place arrows indicating axis
        text(minMC(2),minMC(3),'\leftarrow ','HorizontalAlignment','right','Color','k','FontWeight','bold','FontSize',14);
        text(minMC(2),minMC(3),['   ',num2str(minMC(3))],'HorizontalAlignment','left','Color','k','FontName','Verdana','FontSize',9,'FontAngle','italic','FontWeight','bold');
        text(minMC(2),minMC(5),'\leftarrow ','HorizontalAlignment','right','Color','r','FontWeight','bold','FontSize',14);
        text(minMC(2),minMC(5),['   ',num2str(minMC(5))],'HorizontalAlignment','left','Color','r','FontName','Verdana','FontSize',9,'FontAngle','italic','FontWeight','bold');
        text(minMC(2),minMC(4),'\leftarrow ','HorizontalAlignment','right','Color','b','FontWeight','bold','FontSize',14);        
        text(minMC(2),minMC(4),['   ',num2str(minMC(4))],'HorizontalAlignment','left','Color','b','FontName','Verdana','FontSize',9,'FontAngle','italic','FontWeight','bold');
        %Write values of intercept
        text(minMC(2),1,num2str(minMC(2)),'VerticalAlignment','bottom','HorizontalAlignment','center','Color','b','FontName','Verdana','FontSize',9,'FontAngle','italic','FontWeight','bold');
        
        %return second axes to the top:
        axes(hthr);
        %Legend
        % FDR plot
        if 0
        f=figure;
        set(2,'Color','w');
        plot(VAL(:,1),VAL(:,5),'r','LineWidth',2);hold on
        plot(VAL(:,1),VAL(:,6),'LineWidth',2,'Color',[0.5 0.5 0.5]);
        ax=axis;
        plot([ax(1),ax(2)],[minMC(6),minMC(6)],':','LineWidth',2,'Color',[0.5 0.5 0.5]);
        plot([minMC(1),minMC(1)],[ax(3),ax(4)],':','LineWidth',2,'Color',[0 0.5 0]);
        plot([ax(1),ax(2)],[minMC(5),minMC(5)],':r','LineWidth',2);
        plot(minMC(1),minMC(5),'ko','MarkerFaceColor','r','MarkerSize',5);
        plot(minMC(1),minMC(6),'ko','MarkerFaceColor',[0.5 0.5 0.5],'MarkerSize',5);
        text(ax(1),minMC(5),[num2str(minMC(5)),'  '],'HorizontalAlignment','right','Color','r','FontName','Verdana','FontSize',9,'FontWeight','bold');
        text(ax(2),minMC(6),['  ',num2str(minMC(6))],'HorizontalAlignment','left','Color',[0.5 0.5 0.5],'FontName','Verdana','FontSize',9,'FontWeight','bold');
        text(minMC(1),ax(4),[' ',num2str(minMC(1))],'VerticalAlignment','top','HorizontalAlignment','left','Color',[0 0.5 0],'FontName','Verdana','FontSize',9,'FontWeight','bold');
        axis(ax);
        legend('FDR','FNR',-1)
        grid on
        xlabel('Threshold value','Color',[0 0.5 0])
        title('Positive (FDR) and negative (FNR) discovery rates')
        end

end


% ----- call in functions -----
function gridover() % grid overlay
%axis([0 1 0 1]);
ax=[0 1 0 1];
G=get(gcf);G1=get(G.Children(1));hold on

for i=1:length(G1.XTick) % X grid
    plot([G1.XTick(i),G1.XTick(i)],[ax(3),ax(4)],'k:')
end

for i=1:length(G1.YTick) % X grid
    plot([ax(1),ax(2)],[G1.YTick(i),G1.YTick(i)],'k:')
end
plot([0 1],[0 1],'k:')
hold off